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I. INTRODUCTION 

The phase diagram and critical properties of the d— 

dimensional random fi e ^rfgWiS ijB^I < iA^Tlp^^r7fP f ^ - l' 3een a 
subject of great interestffi|ffl'M'BHO'BBli!. This 
is partly beca use-af-its. i TQle^'l .^|d£scribing the striking 
behaviouimMmeeyaMNe^ of diluted anti- 
ferromagnets in a uniform field. However, the model has 
become celebrated in its own right because of the way it 
captures in an extreme way the effects of randomness and 
frustration; and because of continuiHg-nuestions concern- 
ing its lower critical dimcnsij3ntJl3'Btlll3, the applicability 
of dimensional reductiorjBu'caaBd nature of the tran- 
sition and critical behaviourll^Jtlia. As is w.eLL-known, the 
lower critical dimension di turns out to beEtti .that orig- 
inally suggested by domain wall argumentsElu'El, namely 
di = 2. This implies a marginal -behaviour at dimension 
d = 2, and it has been arguec£SHOE3 that here do- 
main roughening effects are important. This present pa- 
per provides a detailed study based on this point of view. 
We explore the consequences of a roughened domain wall 
picture, and at the same time confirm its validity by ob- 
serving characteristics of the basic ingredients as well as 
consequent predicted behaviour in data from numerical 
studies. 

The domain analysis is most conveniently made in a 
bar geometry. The dependence on a bar width L may 
be used to arrive at phenomenological scaling tranforma- 
tions. These in turn give the phase diagram and critical 
properties. In the particular case of low random field 
strength h and temperature T, the domains span the bar 
width, and are well separated along the bar, so the ba- 
sic "flat wall" description (without roughening) is very 
simple. As expected, it results, for d = 2, in a zero tem- 
perature fixed point at which the field scaling is marginal. 
Adding domain wall roughening breaks this marginality, 
and this ingredient is essential at and near d = 2. Con- 
sequences are: (i) the bulk correlation length behaviour 
£ ~ exp {A/Ki) at T = in d = 2, and (ii) for d = 2 + e 
the phase boundary joining the T ^ 0, h = fixed point 
to one at T = 0, h = h* ~ e (1/7) where £ - \h- h*\- 2/e i. 
7 is one of four exponents occurring in the h and L de- 
pendence of the wall roughening free energy and charac- 
teristic scale, and is predicted to be 7 = 2 in the two 
dimensional lattice. These results have been briefly re- 
ported elsewhereEj. 

This description has been here tested by direct numer- 
ical investigations in the bar geometry for d = 2 (strips). 
This is carried out by transfer matrix techniques (an- 
other reason for using the bar geometry) and by Monte 
Carlo analysis using a new thermalization techniqueo. 
In addition, we use data fawn a max-flow algorithm for 
constructing ground statescj which has been adapted for 
strips. 

The domain description involves free energy T and do- 
main size S as well as roughening characteristics. The 
largest eigenvalue of the transfer matrix gives accurate 



numerical data for and both T and 5 are available 
from the Monte Carlo analysis, for comparison with the 
theory. The wall roughening affects both T and 2, but 
its characteristics are most directly seen in measures of 
the wall profile itself. These and their associated expo- 
nents are best provided by the ground state algorithm. 
The resulting comparisons of theory and numerical data 
give a very complete test of the theoretical description 
and convincing support for its validity. 



The outline of the paper is as follows. In section [I A 
the domain picture is introduced and the flat wall theory 
is provided at T — for the free energy, domain size, and 
correlation length of bars, and via ph enom enological scal- 
ing, for the bulk criticality. Section II B generalises the 
flat wall description to low temperatures, and section II C 
describes the numerical approaches and the comparison 
of free energy and domain size data with flat wall the- 



ory. 



Sccti' 



III A describes the domain wall roughening 



at T = 0, for lattice and continuum models via a simple 
approach and a field theory. First, a single "decoration" 
is discussed, then decorations on all scales, to provide 
the modified free energy and domain size, and hence the 
scaling and criticality. This description is generalised to 
low temperatures in section III B| and compared with the 
numerical domain size and roughening data, includi ng ex- 
ponents, from the ground state algorithm in section III C . 
Section |^ states the main conclusions of the paper. 



II. FLAT WALL 
A. Zero Temperature Theory 

1. Introduction 

In this section, we establish the domain wall descrip- 
tion of the RFIM by taking the zero temperature case 
and assuming flat walls. This is the springboard for all 
the further analytic work. Surprisingly perhaps, the re- 
sults of even this very simply case are very good in certain 
ranges of dimensionality, but they indicate also where a 
more comprehensive picture involving roughening is re- 



quired, and this is provided in section III 



The procedure introduced here is to obtain the dis- 
positions of walls by minimizing an energy which is dis- 
cussed below. The resulting domain size is then used in a 
phenomenological (finite size) scaling scheme to provide 
the RG transformation and hence the critical properties. 
Both these steps are most easily carried out in the bar 
geometry discussed below. 



2. The Domain Picture 



Se suppose, following the prescription of Imry and 
, that the contribution of the field to the energy of 



a domain goes like hv V where V is the number of spins 
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enclosed within a domain and h is the standard deviation 
of the field distribution P, assumed to have zero mean. 
This is plausible for any distribution of fields that tends 
towards a Gaussian distribution in the sense of the cen- 
tral limit theorem. Also, the energy change per domain 
due to the exchange interaction is exactly proportional 
to the perimeter of the domain A; so 



TiT = 0) = n d (2JA - c a hVV) 



(1) 



where n& is the number of domains. In this equation A 
and V represent an average perimeter and volume mea- 
surement for all the domains. As it is a constant, the 
basic ferromagnetic energy ~2JNq (where q is the co- 
ordination number of the lattice) is neglected here and 
throughout this paper. The constant Co is of order unity 
and represents both the statistical fluctuations and the 
selection effects of the domains. We describe the origins 



of this constant, and determine its value, in section II C 5 



We have chosen the strip geometry in which to apply 
the zero temperature free energy (Q) for the following rea- 
sons: (i) convenience for application of finite-size scaling 
procedures, (ii) amenability to domain wall arguments, 
and (iii) we wish to later make contact with transfer ma- 
trix calculations. With these considerations in mind, we 
choose a d dimensional cubic system of finite extent in 
d — 1 directions and infinite in the remaining dimension. 
Our analysis gives a length scale El{Ii) which diverges 
like the correlation length £z,(/i) at the fixed point h*. 
We can then determine the bulk critical properties via 
an RG equation of the form h — > h' = R(h) which arises 
from the phenomenological finite-size scaling ansatz 



q(fl(/Q) = 6lW 

L bL ' 



(2) 



Standard RG procedures provide from Rb{h) the critical 
condition and exponents (from the fixed points and the 
eigenvalues, respectively). 



3. The Zeroth Order Theory 

Given the strip geometry described in the preceding 
section, the leading order modification of the ferromag- 
netic ground state where all the spins are aligned is a 
splitting of the system into rid domains. These domains 
are taken to have flat walls that span the width L of the 
strip and are of typical length El, as shown in Figure ||. 
For this picture to be consistent we need El 3> L, which 
will turn out to hold if h is sufficiently small. In this cal- 
culation, and in the rest of this paper, we set our units 
so that the lattice constant is unity. Using the example 
of equation (Eh we may write the excess energy of this 
domain state as 



T(T = 0) = n d (2 JL d - x - c Q h v /E L L d - 1 



(3) 



where rid — N/El. At zero temperature, this energy is 
equal to the Hclmholtz free energy of the system. So, 
to find the equilibrium free energy we extremalize with 
respect to El yielding the solution 



E, =L 



d-l 



cph 
4J 



(4) 



Applying the phenomenological scaling equation (g) 
gives an RG equation for h of the form 

h' = hir^. (5) 

This implies that there is a fixed point at h — which is 
unstable for d < 2 and stable for d > 2. The canonical 
RG prescription combined with (^|) gives the critical ex- 
ponent v = 2/(2— d) for the unstable fixed point at h = 0. 
This result holds for all d < 2. However, it also warns 
us that d = 2 is a special case where the scaling field 
h is marginal. This is connected with the fact that the 
lower critical dimension of the RFIM is d/ c = 2. Higher 
order terms in the free energy, however, can break this 
marginality. The wall roughening decorations described 
in section III achieve this. 



4- Relation of the Domain Size to Correlation Length 

In the previous section we introduced the average do- 
main size, El. This is related, as we now show, to the 
correlation length £l of the flat wall model. 

Consider the configurationally averaged correlation 
function G(r) between spins in two columns separated 
by r lattice units. Suppose that the flat domain walls 
are independently randomly distributed, so that p is the 
probability that a domain wall lies between two arbitrar- 
ily chosen adjacent columns. Then, since all the spins in 
a given domain are aligned with each other, 



G(r) = (I - p)G(r - 1) - P G(r 
= (1 - 2p)G(r - I) 

and G(0) = 1. It follows that 
G(r) = (I - 2 P y 



1) 



(6) 



(7) 



where the last relation defines the correlation length £l- 
However the definition of p shows that its inverse is the 
average domain size El'- 



P = eS 



(8) 



So we conclude that the correlation length is related 
to the domain size by 



In 



f 



1 



(9) 



In the limit of large El this gives £l = El/2, which is 
the same as provided directly by a continuum version of 
the analysis given above. 
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5. The Correlation Length at T — 

If we assume a flat wall picture, then we can actually 
find the zero temperature correlation length of the RFIM 
on a strip exactly, without employing the more intuitive 
energetic arguments presented above. Us ing th is result, 
we c an ch eck the results of both section [I A3 and sec- 
tion II A 4 . This calculation generalizes the work done by 



Farhi and Guttman for the one dimensional RFIM. 

The key is to consider the connected correlation func- 
tion 



Xll = (o'iCTi+l) — 



(10) 



where we use angled brackets to denote the usual ther- 
modynamic average and cr, and a^i are any two spins in 
columns i and i + I, respectively. We remark that when 
T = the thermodynamic average becomes an average 
over the possible ground states of the system. The trick 
to finding the dominant behaviour of this function arises 
from a desire to devise an algorithm that will solve for 
the ground state of the RFIM. In considering the ground 
state, one notes that even if the field is small, the energy 
to create a domain wall, 2 JL, can be accumulated from 
fluctuations in the random field over large domains. This 
is the reason why there is no long-range order for finite h, 
even at zero temperature, rin two dimensions. Note that 
the work of Imry and MaB demonstrates that in d < 2 
the L dependencies are such that this result persists in 
the L — > oo limit. 

With this in mind, let us define the quantity 



k+l L 



(ii) 



-k i=l 



where we temporarily adopt the notation that i labels 
the vertical coordinate on the strip, and j the horizontal 
one. Then, whenever \S(k, l)\ > AJL, it is favourable to 
insert a domain wall at positions k and k + l. When we 
attempt to employ this algorithm on a variety of field 
configurations, one immediately notices that the ground 
state is still not completely determined. 

A region between two opposite random fields, and with 
endpoints k,j such that S(k,j) = (and j is the site 
closest to k when this is so) is a "floppy domain" (FD) in 
the sense that in an average over the ground states, the 
spins in these domains may point up or down, and are 
not fixed like the spins in other regions ( "rigid domains" 
(RD)). 

One also notices that there is increased degeneracy if 
2JL/h is an integer. This is related to the fact that the 
conditions for both || and || FD can be simultaneously 
satisfied. To avoid the additional difficulties involved in 
the analysis of this case, we shall assume that 2JL/h is 
not an integer. 

It is the degeneracy in the ground state that gives rise 
to non-vanishing connected correlations. Indeed, it is the 
FD's that contribute exclusively to xu- For instance, if 



either spin ex, or spin cr i+ ; is in a RD then (cr,-CTj_|_;) = 
(<7i)(<7i+i). Furthermore, if the two spins are in different 
FD, then they are uncorrelated random variables, and 
their contribution to \u a ls° vanishes. Thus, one need 
only consider spins that are in the same FD. 

Thus, we are left with a combination of three effects 
to consider 

• The probability Wi (/) that there exists a FD of size 
l> I. 

• Given that spin Oi is in a FD, the probability W2 
that the spin ai+i is also in the FD. 

• The thermal average and the average over random 
field configurations subject to the constraint that 
the spins <7i and Oi+i are both in the same FD. We 
label this probability W3. 

We also introduce the variable 



C 



2JL 



where [x\ denotes the integer part of the non-integer x. 

We can view the random field configuration on a strip 
of width L as a random walk in 1 + 1 dimensions. At 
each site along the horizontal direction of the strip, the 
walk changes by a height equal to the columnar sum of 
random field values at that point. For instance, if we 
have a binary distribution of fields, then for L = 3 the 
height changes by ±3, ±1, etc. Using this view, 5(0, 1) 
represents the height of the random walk at x — I (see 
Figure||). 

If the aggregated change in height of the random walk 
ever reaches C then we have accumulated enough fields 
in either the up or down direction to form a domain. If 
there are regions between these peaks where the change 
in height is zero, then these are FD. Figure || should also 
make clear the added ambiguity that arises if 2JL/h is 
an integer. 

This random walk picture is very useful in evaluat- 
ing the three previously introduced probabilities. For 
instance, Wi(l) is the probability that the random walk 
returns to the origin after taking / steps. This is known 
to fall off exponentially with I, 



Wi(/) ~exp 



-I 



(12) 



where £l(£) is a characteristic decay length. We shall 
return to an exact determination of this in a moment. 

From simple geometrical considerations, W2 ~ 
For the third probability, we note that the average over 
ground states corresponds to an average over domain wall 
positions. Since the number of places for a domain wall 
to be inserted between two spins separated by a distance 
I goes like I , we could estimate the conditional probability 
W3 as ~ l/l 2 . These are only rough estimates, the prob- 
abilities depend also upon the ratio h/J. But we only 
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need here the I dependence, and the important feature of 
W2 and W3 is their power law decay in I. We conclude 
that the only exponential contribution is from W\ (I), and 
this gives the exponential decay of x, so then the £l(£) 
defined by ( |l2| ) is actually the correlation length. 

Now W\(l) is the probability that a directed random 
walker in 1 + 1 dimensions returns to the origin after I 
steps, without hitting walls at ±£. We solve this problem 
by reference to the transfer matrix formalism of random 
walkstiH. This involves writing down a transfer matrix 
T(n, m) with n, m G {1, 2, . . . , £} which represents the 
probability that the random walker makes a step from 
position 771 to n. In terms of the RFIM, this means that 



T(n, m) 



Prob(^ hi 



(n — m)h). 



(13) 



i=l 



The leading eigenvalue of the transfer matrix, X(C) is 
related to the survival probability by 

Wi(0~A(£)' (14) 

so that the correlation length of equation (JlJ) is 

-1 



lnA(£)' 



(15) 



For L = 1, the transfer matrix is (using (p_ 3[) ) tridi- 
agonal and symmetric, with zeros on the diagonal and 
1/2 on the neighbouring off diagonals. We can generalize 
this result for all L, by noting a general pattern in a row 
through the middle of the transfer matrix: 



(i-j) 


...-4 


-3 


-2 


-1 





+ 1 


+2 


+3 


+4... 
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16 


16 


16 


16 ■ ■ ■ 



The numerators of the elements follow the same pat- 
tern as the numbers in Pascal's triangle, and the denomi- 
nators are simply 2 L . This is due to the Gaussian nature 
of the random walk whereby T^,(n, m) = \Ti(n,m)] L . 
These matrices are all diagonalized by the same type of 
vectors: 



4\ 



sin (jirq) . 



(16) 



where the wavenumber q is chosen to be consistent with 
the boundary conditions of the matrix, i.e. q = n/(C+l). 
The pattern of entries in the transfer matrix gives the 
eigenvalues of these vectors as 



A = cos L (<77r). 



(17) 



Clearly, the largest eigenvalue is the one with q = q\ = 
1/(C + 1). If h is small, then L is large and q\ is small 
so that equation (|l^) gives the correlation length as 



L In cos 
2J 



2L 



7l7> 



2J/h 
2 



2 , 

16tt 2 



•0(1). 



(18) 



It is easy to generalize these arguments to higher dimen- 
sions via the replacement L — ► L d ~ x . Then ( |l8|) agrees 
to highest order in h with our earlier result for the T = 
small h correlation length (from combining (|J) with (|J)). 
However, it is important to note that this result came 
from an analysis of the connected correlation function, 
whereas our prior analysis applied to the disconnected 
correlation function. Although the correlation lengths 
derived from these two quantities should diverge in the 
same manner, they may differ by a multiplicative con- 
stant ( e.g. Cq). Indeed, this turns out to be the case in 
section 



IIC 



Although this calculation also assumes a flat wall sce- 
nario, it should be accurate when £^ 3> L and is exact 
for L = 1. 



6. Conclusions 

In this section, finite-size scaling in the strip geometry 
has been used to study the zero temperature fixed point 
of the RFIM. We found the critical exponent v = 2/(2 — 
d) for d < 2. This result indicates the need for a more 
sophisticated analysis when d = 2, and this follows in the 
subsequent chapters. 

The procedure involved the determination of the free 
energy and the typical domain size and its relationship 
to the correlation length. These results are confirmed by 
a direct calculation of the correlation length of the RFIM 
on a strip, assuming flat domain walls, which agrees per- 
fectly up to the coefficient cq. 

This provides a consistent treatment of the RFIM at 
zero temperature and low random field strength away 
from the lower critical dimension. However, when d = 2 
fluctuations are very important to the critical behaviour. 
So, the next section studies the effects of fluctuations in 
the domain wall shape, and this will complete the de- 
scription of the zero temperature RFIM near the critical 
point. Comparisons with numerical work and investiga- 
tion of thermal effects will be provided in later sections. 



B. Low Temperature Theory 

The extension of the domain wall analysis to non-zero 
temperature requires the construction of an entropy. The 
procedure will then be to minimize the free energy T = 
U — TS, and so determine the characteristic lengths for 
exploitation via phenomenological scaling. By working in 
the low field and low temperature limit, i.e. low H and 
T, we can use the well separated strip-spanning domain 
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pictures developed in the previous section. For the zeroth 
order (flat wall) description, the entropy then turns out 
to be quite trivial, corresponding to the number of ways 
of laying down flat domain walls with an average spacing 
defined by Ex = N/na- Thus, the entropy is 



So (Ex) = In 



{n dJ 


= N 







ln(S L - 1) 



-In (1-3 



(19) 



choosing units so that ks = 1 and using Stirling's ap- 
proximation in the thermodynamic limit. The free en- 
ergy generalizing ([!]) to T ^ is then 



T = n d { 2JL d 



Extremalizing with respect to 
equation for the domain size 



WSx^- 1 ) -TSo(Sx). (20) 
Ex gives the following 



2JL d ~ 



C h T ±zl„l/2 

2 L 



Tln(Sx 

The various limits of this equation are 



1) = 0. (21) 



Ex = 1 + cxp (*y£-^ , h = 0. 



(22) 



As a check on equation (|22j), we can use it to investigate 
the low temperature behaviour of the pure Ising model 
(zero field) . By using this form for Ex in the phenomeno- 
logical renormalization equation (^) we get the following 
RG equation for T j J 



W- 1 (t) 



2(bL)*-l-Qnb){$) 



(23) 



This has a fixed point at T — which is unstable for 
d < 1; it is stable for d > 1, implying the existence of a 
finite critical temperature which we know the pure Ising 
model should have in d > 1. In d = 1, ( f22| ) gives the 
usual pure Ising result £ ~ exp(2J/T). 

Thus, in d = 1 the two scaling variables (scaling like 
a length) are exp(2J/T) and (from (§) h- 2 ^ 2 - d \ It 
follows that in d = 1 the correlation length can be written 
in the crossover form 



£ = e 2J ' T 4>{h 



2 e 2J/T 



(24) 



where 4>{x) — ► c a or c&/x for a; — > or oo, respectively, 
with c a ,Q, constants. 

For g? > 1 the stability with respect to tempera- 
ture of the (T, h) = (0, 0) fixed point means that no 
crossover form then applies at low T, h, but instead 
j l 2/(2-d)rpi/(d-i) j g mv ariant under scaling. 

Because of the attractive thermal nature of this T = 
fixed point at small h and d > 1, the thermal scaling has 
little influence on the low T low h critical behaviour for 



d > 1. Nevertheless, thermal effects contained in ( |20| ) 
and ( ^l| ) can substantially modify the behaviour of the 
free energy and domain size in the low T, low h regime 
for d > 1, masking the field dependences more significant 
for the RG scaling in this regime. So ( pfj| ) and ( |2l| ) will 
be required for interpretations in the next section, where 
numerical techniques are used to confirm the basic scaling 
picture so far discussed. 



C. Numerical Evaluation 

1. Introduction 

We wish to test the results developed in the previous 
section by comparing them to independent numerical in- 
vestigations. If we find reasonable agreement between 
the two, then one can feel confident about the veracity 
of both results. 

The structure of the analytical approach is such that 
the numerical tests can be applied at various stages. The 
most obvious test of such a theory is on its final scal- 
ing predictions for the bulk system criticality. However, 
this is difficult since it requires the investigation of a sys- 
tem large enough to show the bulk critical behaviour. It 
is much easier to apply the numerical techniques to the 
finite size strip geometry, as is done here, and to investi- 
gate the non-critical ingredients from which the scaling 
transformation is built, and hence establish the applica- 
bility of the basic flat wall procedure. 

Three techniques have been found to be successful, and 
their results are reported in this paper. There are a trans- 
fer matrix calculation, a modified Monte Carlo simula- 
tion, and a ground state algorithm. 

The transfer matrix method is appropriate to the strip 
geometry, and can in principle provide the phenomeno- 
logical RG transformations directly and very precisely, 
as has been exemplified for many low dimensional lat- 
tice systems, particularly those with non-random trans- 
fer matrices. 

Very accurate results have also been obtained for ran- 
dom bond and site diluted Ising models and a modifica- 
tion of those techniques is used hereE3. Because of the 
quenched disorder present in such models (and in the 
RFIM), the first (dominant) Lyapunov exponent of the 
transfer matrix product gives the average free energy of 
the strip, which is one of the two key ingredients in our 
domain wall analysis which we would like to check by the 
transfer matrix and other techniques. The other ingredi- 
ent is th domain wall size or correlation length. Typically, 
the correlation length is provided by the TM approach 
via 
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-l 



ln(A!/A ) 



(25) 



where Ao and Ai are the largest and second largest eigen- 
values, respectively. However, this relation gives the most 



G 



probable correlation decayed, which need not be the same 
as the average correlation decay. In fact, the correlations 
may follow quite complex distribution functional!, and 
the difference between "most probable" and "average" 
decays becomes especially important when the effective 
interactions may diffe r in s ign, such as in spin glasses or 
random field problemscBcIJ. For this reason, the compar- 
ison with transfer matrix data will be restricted to the 
free energy. 

Monte Carlo simulation is the second technique used. 
The usual equilibriation problems set by critical slowing 
down are here compounded by the effects of very large 
energy barriers in situations (d < 2) where the criticality 
is at T = 0. The problems are partly alleviated by the 
strip geometry, particularly in d — 2, where character- 
istic lengths of order h~ 2 (see (Q)) replace exponentially 
large correlation lengths of the bulk system (see (0)). 
However the equilibriation times remain excessively, long, 
and a new algorithm (described in more detail intfl) was 
devised and employed to markedly improve thcrmaliza- 
tion times. Briefly, this algorithm continues the ideas of 
multigrid methods with self adjusting renormalisations of 
parameters to allow for changes to be made on all length 
scales via moves of blocked spins. This permits the calcu- 
lation of average energies, domain sizes, and correlation 
lengths, for comparison with the analytic predictions. 

The 6-t\al numerical technique used was an exact al- 
gorithmcd to construct the ground states for both bi- 
nary and Gaussian random field configurations. Ground 
states for models with arbitrary random fields and ar- 
bitrary but not frustrated exchange interactions may be 
found in polynomial time by mapping the optimization 
problem to_a min-weighted-cut problem on an associ- 
ated graphcj. A new min-cut-jnax-flow algorithm!!!! was 
first implemented by OgielskS to demonstrate the prac- 
ticability of this method. We here apply the algorithm 
to the strip geometry, and hereafter refer to this method 
as the "max-flow algorithm". From these data it was 
possible to obtain T = properties including average 
energy, domain size, as well as domain wall roughening 
characteristics which will be discussed later. 

We now proceed to the numerical data which test the 
flat wall picture discussed so far. The transfer matrix 
determinations of the free energy are first described and 
compared with the analytic work, followed by the Monte 
Carlo results and comparisons, then max-flow results and 
comparisons. 



2. Transfer Matrix 

The simulations were done on strips of width L = 
2, ... ,7 and length N — 10 5 and a binary random field 
distribution. Periodic boundary conditions were applied 
in the "finite" direction whereas all possible trial state 
vectors were applied to termini of the strip at each end 



of the "infinite" direction^ This procedure was then re- 
peated for three different realizations of the random field 
configuration, and the resulting free energy 



F/N + 2 = -fcTlnAo 



(26) 



averaged over these realizations and over the 2 L trial vec- 
tors. The average was calculated in the root-mean-square 
manner, and the associated error in the free energy was 
calculated as the rms error in the free energy at each 
random field configuration. 

We wish to test the predictions of the zeroth order flat 



wall theory (section II A), in which the free energy at 
very low temperature and d = 2 is 



T Q = — 



Nc 2 h 2 
8LJ 



(27) 



To this end, we have generated data for L E (2, 7) at 
T = 0.1J. For lower temperatures, the elements of the 
transfer matrix become larger than the computer can 
handle. If the strip width is much larger than 7 the run 
times become unreasonably long (over a week). Never- 
theless, these strip widths are adequate for our purposes. 

As a test of this functional form, we plot \ii(!Ftm / J) 
versus ln(/i/J) for various strip widths L, where Ttm 
represents the free energy determined by the transfer ma- 
trix method. Straight lines of slope 2 are expected, and 
this is indeed the demonstrated in Figure EjL 

Straight line fits to these plots yield the slopes given 
in Table J. We include the value of the Q statistic, which 
is a measure of the probability that a value of x 2 as poor 
as the value we have found should occur by chance. Q is 
always in the range CO. 1) and a value of Q > 0.1 usually 
denotes a reliable fitEl 

The agreement with the predicted slope of 2 is ac- 
ceptable, with increasing accuracy at larger L. How- 
ever, there appear to be systematic errors, as one can see 
from a close inspection of Figure |3[ This is an indica- 
tion of the importance of higher order decoration effects 
and finite-temperature effects, occurring even at these 
low temperatures and high fields. 

We now attempt to obtain a best value for the constant 
Co, so that we can compare the predicted free energy with 
the numerics, and extract subdominant features noted 
above. 

A simple fit of the numerical free energy to equa- 
tion (^tJ) is statistically unreliable, and gives residual L 
dependence in Co, and other artifacts. We concluded that 
this was due to the data being at T = 0.1 J and not zero 
temperature. Indeed, for small h, the leading temper- 
ature dependence can be very significant. The zeroth 
order finite temperature theory of ( |20| ) and (^l|) implies 
that for large Si in a two-dimensional lattice 



As with all the numerical data presented in this paper, these 
calculations were performed on DEC Alpha workstations. 
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T = -N- 



8LJ 




L 



4J 
c h 



(28) 



Since the logarithm diverges for small h, the tempera- 
ture dependence becomes increasingly important in this 
regime. 

To fit the free energy data properly, we need a free en- 
ergy form applying for all h, T as long as Si is large. This 
can be obtained by reparameterizing the relations ( pOj ) 
and ( |2l| ) giving the finite temperature flat wall theory, 
provided the theory is correct for large Si,, as we have 
argued. Using h, T small, large, d = 2, the approxi- 
mate forms to rcparameterise are the excess free energy 



T 
~NL 



f 



-c h 



2VSIT 



and the extremalizing equation 



,hL 1/2 
IT 



-1/2 



1 



In S 



il exp 



-2JL 
T 



(29) 



0. (30) 



We now define a new variable y — c(T) V^l, where 
c(T) = exp(— LJ/T). Then, the intensive free energy 
is 



f 



(31) 



8(3yT 

and the minimizing equation for S^, has the simple form 
lny = -A/ (32) 

where 

cohVX 



AcT 



(33) 



We may then eliminate y between ( J3l| ) and (g2|) to 
arrive at an equation for / involving (3. Inserting its 
explicit form leads to 



LJ = f -c fe \ c 2 h 2 
T U \2fVLj 8/T' 



(34) 



This is an implicit equation for / as a function of h, T, L. 
We will use it to solve for / numerically and so determine 
cq. The usual x 2 statistic was calculated 



x v v * J 



(35) 



where Ci is the estimated error on Jztm- This statistic 
was minimized by a bisection routineo, which also calcu- 
lates the Q statistic. 



As well as obtaining Co, we can test whether non- 
Gaussian effects occur for small L, or perhaps the flat 
wall description breaks down for large L by varying the 
range of the fit. The values of T used were 0.1J, 0.2 J 
and 0.3 J, with h £ {0.05,0.25}. Over this range of T 
and h, the zeroth order theory predicts large correlation 
lengths at least > lOi for even a liberal estimate of Co, 
which should make (|34| ) applicable. This assumption was 
also demonstrated to be self-consistently true by using 
the fitted value of the selection constant to predict Hj> 
Over this range of fields, then, the results are presented 
in Table ||- 

One notices the striking improvement to the fits when 
the lower L values are dropped (Q = 0.686 is remarkably 
good). We attribute this to non-Gaussian selection ef- 
fects at lower L values, where the domains do not sample 
over as many field configurations as a domain of the same 
size at large L. Since the number of field configurations 
for a column of length L is 2 L , we would expect a marked 
decrease in this effect as we increase L. 

The value c ~ 1.85, found by fits of the theory to 
the free energy data, does n ot agree with that (co = 
7r) predicted in section II A 5 by a random walk analysis 
of the disconnected correlation function. As explained 
this is due to the different correlation 



IIA5 



m section 

functions used in obtaining this (non-universal) constant. 

The large Q statistic gives much confidence in the form 
of T predicted from the flat-wall ansatz. 

Figure [I] illustrates the quality of the agreement be- 
tween the predicted form of the free energss and that pre- 
dicted by equation ( |34"| ) using Co = 1.85.B. One should 
remark that, while the fits were done for h S (0.05, 0.25), 
the predicted free energy seems to lie within error for 
a much larger regime, especially for larger strip widths. 
Indeed, for L = 8,9 the predicted free energy lies within 
error of the transfer matrix data all the way out to 
h = 1.0 J. 

Transfer matrix data were also taken for h/J £ (1,6) 
and low T, to determine the effect of this increasing field 
strength. Typical results are shown in Figure | One no- 
tices that there are three separate regimes in the graphs: 

• h/J £ (0,~ 2). In this regime, the free energy 
seems to go roughly like the square of the field. As 
emphasized above, this dependence is close to that 
predicted by the flat wall theory, with increasing 
accuracy as the strip width, L, increases, which 
is important for our scaling discussion of critical 
effects. 

• h/J £ (« 2,4). Here, the free energy crosses over 
from quadratic to linear dependence on h, due to a 
rapidly decreasing domain size. 



The more efficient Newton-Raphson type routine could only 3 Th H t „ t T _ n 2 r whilp it was „ spd - m fhp fits disri i av5 , 

e used it we had numerical derivatives of f. these are not c , J 

., r i „ ,. XTl ■' no new features, 

easily found from equation (pj). 
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• h/J > 4. The field strength is greater than the 
energy cost of flipping a spin, and so all the spins 
align with their local field. Thus the free energy is 
linear in h, with unit slope. 

To sum up: the transfer matrix approach has pro vided 
a ri gorou s test of the energetic arguments of section [I A 3 
and II B . It provided excellent agreement over the range 
tested, giving confidence in the applicability of the the- 
ory. 

However, a comparison of quantities like correlation 
length, domain size, etc., would be a more rigorous check 
and would also provide a better picture of the physics 
underpinning the results. Secure results on these lengths 
are not so far available from the transfer matrix method 
and so we turn to the Monte Carlo analysis, which can 
provide them. 



3. Monte Carlo 

Using the Block Monte Carlo algorithm described in 
referencetJ we now investigate other physical quantities 
which are directly acessible via our theoretical frame- 
work, but not easily obtained by transfer matrix analysis. 

a. The Average Energy: First, we compare the con- 
figurationally averaged energy found by the Monte Carlo 
rou tine to that predicted by the flat wall results of chap- 
This quantity is the first term in (pol) and fri- 



ll E 



ter 

volves Sx which is given by numerical solution of (21). 
The theoretical and Monte Carlo data points are shown 
for L = 1 and 4 in Figure ^. 

Similar results are obtained from the data for L = 
2. Fitting the data to the theoretical form (|3|) gives 
Co = 1.75 ±0.05, in reasonable agreement with Co = 1.85 
need ed to fit the transfer matrix free energy data (sec- 
tion 



IIC2) 



The theory agrees quite well with the Monte Carlo 
data up to a temperature, that diminishes with increasing 
field strength, at which Sx ~ L and the criterion for 
the validity of the theory breaks down. Yet, there is 
agreement within error up to T = 0.7 J even for h = 0.5 J 
when L = 4. This is a further indication that the zeroth 
order theory captures the essential flavour of the RFIM. 

b. The Domain Distribution and Correlation Length- 
Domain Size Relationship: A key assumption in our 
derivation of the zeroth order entropy (section [I B ) and 



a t ion 



of a relation bet ween d omain size Sx and corre . 
length £x (section II A 4 ) was that the domain walls were 
randomly placed on the strip. This can be tested ei- 
ther by investigating the domain size distribution, which 
should then be Poisson, or by exploring the predicted 
relation between Sx and £x- 

Since the standard deviation and mean of a Poisson 
distribution are the same, we have evaluated these two 
quantities for L = 1, 2 and 4 for h < 0.5 J and T < J. 
In every case they are found to be the same within error 



(never more than 5%). While this result is not definitive, 
it is very suggestive of a Poisson distribution. 

Equation (||) expresses the relation between the cor- 
relation length £x> and the domain size Sx when on a 
lattice, for randomly distributed walls. Figure fj] shows 
a test of this relationship for L = 2. Very similar plots 
arise for L = 1 and 4. In each case, the agreement with 
the predicted form is remarkable, and gives us confidence 
that, at least fo r larg e correlation lengths, the methods 



used in section II A 4 are correct and the domain walls 



really are Poisson distributed. 

The final support for this conclusion comes from direct 
analysis of the probability distribution of domain sizes it- 
self for L = 1, 2,4; h = 0, 0.2 J and T up to J. The results 
for L = 1 are shown in Figure The resulting size distri- 
bution is very close to exponential (providing strong evi- 
dence that the domains really are Poisson distributed) for 
h = 0, at least for smaller L's and larger temperatures. 
However, when h = 0.2J and T = 0.4J the distribution 
seems to decrease at lower S. This is evidence for do- 
main wall repulsion, which is expected at higher fields 
when the domain walls begin to approach and undergo 
appreciable roughening (see section III ) . For h < J , how- 
ever, this effect is only noticeable over about the lower 
2% of the domain size distribution. Otherwise, the dis- 
tribution is still exponential. So, the Poisson assumption 
appears to be approximately correct, and gives the cor- 
rect behaviour at large domain sizes. Thus, the flat wall 
entropy should still be a valid starting point, although 
a more complete description is obviously required for a 
full understanding, particularly in the large size regime 
of interest for scaling. 

c. The Domain Size Itself: Having justified assump- 
tions implicit in the flat wall theory, we wish to test its 
predictions further. In this section, the domain size Hx, 
predicted by (|2l]) is compared with that measured by the 
MC routine. 

The results are shown for L = 1 and 4 in Figure |^. 
The data for L = 2 looks similar. The flattening of the 
data at low T's is due the finite length of the simulated 
system; data in the flattening regime are discarded in 
quantitative comparisons. The constant cq was deter- 
mined by a best fit to the data at T = 0.1J, where we 
might expect the flat wall theory to be accurate. The re- 
sult, Co = 1, contrasts with Co = 1.75 ± 0.05 given by the 
average energy fits presumably because selection effects 
act differently for length scales and energetic quantities. 

Contrary to the previous convincing agreements with 
theory, the comparison of the predictions of ( |2l| ) to the 
domain size data is relatively poor. The reason for this 
poor agreement is not well understood, but may be due 
to non-Gaussian effects at smaller strip widths. How- 
ever, we do get reasonable agreement (discarding flat- 
tened data) at the lowest of temperatures, and for larger 
L. These are the conditions that matter most for the 
scaling constructions, where the correlation length and 
domain size at T = and large L is of particular im- 
portance. These are discussed below using data from the 







ground state algorithm. 



4- Max-Flow algorithm 

We have used the max-flow algorithm^! to find ground 
states and hence ground state properties for many ran- 
dom field configurations generated from binary or Gaus- 
sian distributions. In the large scale limit, these two 
distributions are expected to give the same results, but 
there are differences, discussed in section [II C| , for small 
system or domain sizes. These differences were not as 
pronounced in the finite T data obtained from the trans- 
fer matrix and Monte Carlo routines described above. In 
the present case, it becomes advantageous to use results 
from the Gaussian distribution as we now do. 

These data were generated on strips of width L = 
2, . . . , 11 and length N = 10 3 . These results were av- 
eraged over 100 independent random field configurations 
and the error bars shown represent the statistical varia- 
tions in the data when averaged over these runs (as op- 
posed to the width of the distribution of a variable quan- 
tity like the domain size). As the computational time 
required to generate a ground state via this algorithm 
grows like n 3 , and the run time for an n = 11 x 10 3 
system was on order 10 days, it was not possible to sig- 
nificantly increase the system size beyond this limit. 

The typical domain size is shown in Figure [l(J Also 
included in the Figure is a comparison with the theoret- 
ical obtained from the random walk analysis (sec- 
tion pi A 5 ) of the correlation length £j, using a simple 
fitted constant of proportionality to convert this to a do- 
main size. A quantitative assessment of the quality of 
the agreement is provided by the correlation coefficient 
r = 0.9968 and the statistic x 2 = 0.426. The fit gives a 
value of Co = 1.98 ± 0.02, in reasonable agreement with 
the values obtained by both the transfer matrix and block 
Monte Carlo results at finite T. 

Also considered was a fit to just the h~ 2 leading term 
in the small h expansion of the theoretical 2/,. This 
corresponds to (Q) (i.e. to the basic T = theory and 
differs from the full by 0(l/h), significant except in 
the low h limit. The comparison of fits to this reduced 
form with those to the full show that if fields up to 
h = 2.0 J are included , the fit to the reduced form gives 
a x 2 — 4.215. This shows a much poorer quality of fit 
than that to the full theoretical form. The neglect of 
the 0(1/ h) terms is one of the sources of the reduced 
quality fit observed in the basic theor y ana lysis of the 
Monte Carlo domain size data (section II C 3 ). However, 
it is difficult to incorporate these higher order terms in 
a consistent finite temperature theory. This shortcoming 
disappears in the small h regime important for scaling. 

This completes the comparisons of numerical analyses 
of (transfer matrix, block Monte Carlo, ground state al- 
gorithm) with the basic theory. The conclusion is that 
this theory provides an excellent account of the free en- 



ergy and average energy; and an account of the domain 
size and correlation length which differs from the nu- 
merical data by effects whose origin is becoming clear, 
particularly so at low temperature. These effects are 
the domain wall repulsion due to non flat wall effects 
(wall roughening), seen in the domain distribution of 
section [I C 3 , and the (understood) 0(1/ h) and other 
terms by which the basic theory differs from the ran- 
dom walk one. For these reasons, and for those given 
in section II A 6 it is necessary to extend the theory by 



including wall roughening effects and to verify their im- 
portance, which is done in the next section. 



III. DOMAIN WALL ROUGHENING 
A. Theory 

1. Introduction 

The effect of domain wall roughening on the free en- 
ergy, and hence the correlation length is estimated here 
by a decoration method applicable in general dimension- 
ality. Lattice and continuum cases will be studied, and 
these turn out to have different roughening characteris- 
tics. We first outline a simple yet effective calculation 
(given previously inES), first on a lattice and then in the 
continuum. Then, a functional field theory in the contin- 
uum is devised which provides a quantitative description 
of the effects of wall roughening. The varying results of 
these three approaches will then be compared. These 
results are next used to perform wall roughening on all 
length scales. This provides an RG transformation in- 
cluding the terms which break the marginality of d = 2. 
Hence, we can derive the scaling behaviour in d = 2, and 
the critical point and critical exponents of the RFIM in 
d = 2 + e. All these calculations are carried out at zero 
temperature for simplicity, and will be later generalized 
to finite temperature. 



2. The Lattice Theory 

On a lattice, the basic perturbation of the domain wall 
("decoration"), to be later introduced on all scales, can 
be taken to be a rectangular shape (see Figure |ll|), of 
width a and height b. The associated excess volume 8V 
and the excess area 6 A are, in d dimensions 



SV = ba d ~\ 8 A = 2(d - l)ba d ~ 



(36) 



We assume that the domains are sufficiently large and the 
domain wall fluctuations sufficiently small (i.e. ^> b) 
that the fluctuations in the wall may be considered as 
uncorrelated from those in the body of the domain. Thus, 
we may treat the decoration as a separate entity from 
the bulk of the domain and apply the free energy (jl]) 
independently to it. 
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Substituting the values ( pq ) into equation ([j]) gives the 
zero temperature excess free energy 

8T(T = 0) = AJ(d - l)ba d - 2 - c a hb 1/2 a {d ~ 1)/2 . (37) 

We minimize this free energy with respect to b to arrive 
at the relation 



b = 



( coh 



,3-d 



\SJ(d- 1) 
This relationship is of the generic form 
b = c(h/J) K a c 



(38) 



(39) 



which defines a wandering exponent £, an d a further 
"field" exponent k. We see that C = 3 — d,K = 2in 
this simple description 

Inserting ( pq ) into ( |37| ) gives the equilibrium value of 
the excess free energy and in particular its power law 
dependence on field and on the scale a, 



For self-consistency, we must require that (0) holds 
at large scales a only for d > 2 and that hj J is not much 
larger than unity, or that d = 2 and ft/ J <C 1. Similarly, 
equation ( [431) only holds if d < 2. We now check the 
assumption that the wall decorations are decoupled from 
the bulk of the domains, i.e. that b(L) <C S^. Using the 
form of Si given by the zeroth order decorations (^) and 
the solutions for b(L) given in equations ( f42| ) and ( |i^ ) 
above we conclude that h/J must satisfy 



f d — 2 

ft- I L~ , if 6 < a 
J I L^r , if b > a. 



(44) 



In either case then, for large L and d > 2 there is a gener- 
ous regime for h/J which satisfies the assumption. When 
d < 2, if L is finite we can always find h/J sufficiently 
small that (fjij ) is satisfied, but the allowed regime for 
h/J decreases as L increases. 

Thus, these simple arguments give the wandering ex- 
ponent of the domain wall, f , 



,3-c 



J 16(d-l) 



(40) 



3. The Simple Ansatz for the Continuum 

We now perform a similar calculation to that in the 
above section but in a continuum rather than a lattice. 
It turns out that these two cases have different exponents. 

By analogy to the previous section, the basic decora- 
tion of the domain wall is described by the typical height 
b and the typical width a (see Figure [y]). These will 
again be related in the manner (|39|). Assuming axial 
symmetry, we denote the excursion of the domain wall 
from its mean position as z = bf(r/a) where /(0) = 1. 

?y (0) is 



Then, the excess surface area in the free energj 



SS(a,b) = /'drr'-* ( yfl + (b&f ^))' 



1 



b 2 a d-3 

ba d - 2 



if b < a 
if b > a 



Constants of proportionality have been disregarded here, 
and will be elsewhere in this section, since they do not 
affect the exponents in which we are interested. 
Similarly, the excess volume is 



SV{a,b) = / drr d - 2 bf[- 



ba" 



(41) 



Using the above forms in the free energy (|]J) and min- 
imizing with respect to b yields the following relations 



ft) 



2/3 5=d 

a 3 , 



-,3-d 



if b < a 
if b > a. 



(42) 
(43) 



c 



3-d 



ifd> 2 

if d < 2 and h < 1 



and the free energy 
5T(a) - ■ 



,4/3 d+1 



if d> 2 

if d < 2 and ft < 1. 



(45) 



(46) 



4. 27ie Fidd Theory 

A more complete treatment of this problem allows the 
entire shape of the interface to be determined by the sys- 
tem, and not just the height to width ratio as in the 
previous treatment. This will provide a check on the 
results of the previous section, and will provide more de- 
tailed information on the average equilibrium shape of 
the domain walls, and on proportionality constants. 

In order to achieve this, we will establish a functional 
equation for the shape of the interface by rewriting the 
free energy of the domain wall as a functional of the 
wall profile P. The equilibrium wall profile is then found 
by extremalizing this free energy via the Euler-Lagrange 
equation. The results of this section only hold in d > 2. 

We use ad-1 dimensional cartesian coordinate system 
placed on the hyperplane where the flat interface would 
have been. Then, the deviation of the domain wall from 
this hyperplane is given by the function P(x). Using the 
model of equation ([!]) , this decoration has a free energy 
at T = of 

F[P(x),VP(x)] =2J yfl + (VPyd^x 



coh J / P{x)d d ~ 



(47) 
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where both integrals are over the hyperplane defined by 
Xi G [0, L]. This free energy is minimized via the gener- 
alized Euler-Lagrange equation: 



d-1 

E 



ST 



ST 



dx t 5{dP(x)/d Xl ) SP{x) ' 

Note that the RHS is a constant, and define it as 

-c ft- 



N J 



4^ '/ P{x)d d - l x 



(48) 



(49) 



This leads to a rather complicated d-dimensional dif- 
ferential equation. 

d-i {afp \i + (vp) 2 ] - EiZldiPdjPduP} 



i=i 



[1 + (VP) 



21 3/ 2 



(50) 



To simplify things we consider the isotropic solu- 
tion, which only depends on the radial coordinate r 2 = 

£•=!*!• 

Then, with d r P(x) = P'(r) the isotropic equation be- 
comes 



Nn = 



(i=2) p'( r ) (1 + P'(r) 2 ) + P"(r) 



(1 + P'(r) 



2^/2 



(51) 



Note that this equation may alternatively be obtained 
by extremalizing an isotropic free energy functional, 



5T[P(r),P'(r)]=2jn / a/1 + P'(r) 2 r d - 2 dr 



c hJn / P(r)r d ~ 2 dr. 



Here 



n 



(d-i) 



(52) 



(53) 



is the area of the d— 1 dimensional unit ball. However, in 
the more general case the isotropic solution is a member 
of the larger set of optimal solutions. 

It may be checked by direct substitution that the so- 
lution to the wall decoration problem encapsulated in 
equations (49) and ( |5l| ) is a hypersphere of radius R 



P(r) = ^R 2 - r 2 



Pn 



(54) 



where R satisfies the requirement —1/R = No. The con- 
stant Pq is determined by the choice of the boundary 
conditions P(r — L/2) = 



Po 



] R 2 



(55) 



In addition, the constraint 
consistency relation 



imposes the self- 



R 2 



J 2 (d-1) 2 A(5V) 



c 2 h 2 



(56) 



which we use to determine R(h). Here SV is the extra vol- 
ume added by the decoration 5V = Q J^ 2 P(r)r d ~ 2 dr. 



Thus, ( |56| ) is difficult to solve for R(h). However, we 
are interested in the regime h/ J < 1 where we expect the 
decorations to be small (R 3> L) and we may expand in 
powers of L/R. To first order in 1/P the excess volume 
and area are 



sv = 



n 



R{d 2 



8A 



2P 2 (d+l) \2 



d+l 



d+1 



(57a) 
(57b) 



Then, equation (pq) gives the equilibrium radius as 



p ,'c h\ 3 ( 4 Q(d-l) 



(58) 



Inserting the first order forms (|57]) evaluated at the 
equilibrium radius (|58| ) into (^) then gives the first order 
free energy fluctuation 

T wall (L) = -c (^y\^+0(h 8 / 3 ). (59) 



Here c is the numerical constant 



fiV3 



(d+l) 1 /3(d_ 1)2/3 



_ 2 -!/3 + 2 8/3 



2 (d+7)/3 



(60) 



and is positive for all d > 1. This shows that wall rough- 
ening is always favourable for d > 2. 
Finally, the wall wandering height is 



P(r = 0)~h 2 / 3 L^+O(h 8 / 3 ). 



(61) 



We now compare these three approaches to the wall 
decoration problem on the basis of their predicted be- 
haviours of the wall wandering height b and the free en- 
ergy of the wall ST '. We define the wandering exponent 
C, the free energy fluctuation exponent u>, and associated 
exponents 7 and re by the following equations 



ST 
b 



h K a c , 



(62a) 
(62b) 



where a is the base length (L in the field theoretic de- 
scription). One sees that all the exponents of the simple 
continuum ansatz and the full field theory agree for d > 2. 
The complete comparison between all three treatments is 
given in Table IV, and the results are shown graphically 
in Figure n3. 
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The predicted exponents in the continuum model agree 
with those of other extensions of the Imry-Ma type ar- 
gument done by NattermantL3. On the other hapd, the 
lattice exponents agree with the work of BindeiEJ, also 
for the lattice. GrinsteinE-3, finds the same values for 7, 
k, and £ as we do, but finds 8 = 1. It is also worth 
noting that the replica symmetric approach of Parisi and 
SourlascEl for the continuous model gave C = (5 — d)/2, 
but the inclusion of replica symmetry breaking gives 
C = (5 — d)/3oc3. Finally, the lattice result is consistent 
with the results of transfer matrix-calculations done in 
two dimensions by Fernandez et aZJ12l and with numerical 



ground state calculations reported in section [II C 



scale by a factor of n l so there are n' l( - d ~ 1 ' decorations 
at this stage, and each gives a free energy contribution 
5J-(L/n l ). Hence, the total free energy of all the decora- 
tions will be 



i=0 v 7 



(66) 



with k as above. 

In this way, then, we can estimate the effect of wall 
roughening on all length scales, and remove the restric- 
tion that the wall shape must be smooth. Such general- 
izations are crucial to get the correct marginality break- 
ing free energy, as we shall see. 



5. Decorations upon Decorations 

The decorat ion to the flat domain walls described in 
section [II A 4 is (to lowest order in hj J), the most favor- 
able twice differentiable form of the domain walls. How- 
ever, discontinuities may be a feature of the "true" solu- 
tion. We may incorporate such discontinuities by allow- 
ing the decoration to be repealed on all length scales, fol- 
lowing a procedure of BinderEI Basically, having made a 
basic change in the shape of the domain wall on scales of 
length L, one then looks at the domain wall on a length 
scale L/n where n > 1 is some arbitrary integer, whose 
choice, we hope, should not strongly affect the results. If 
the random field is small, the wall will be slowly varying 
on such a scale and may be viewed as a hyperplane, or to 
preserve the simplicity of the discussion by taking d = 2, 
as a straight line. This straight line is then decorated 
in the same manner as the original wall. This process is 
repeated on length scales of order <n = L/n l , i = 2, 3, . . . 
until a cutoff length scale is reached (see Figure |l3|) . Re- 
verting to the case of general d, this cutoff occurs at a 
stage k where the smallest length scale of decorations is 
the lattice constant (unity), i.e. 



min (a k ,b k ) = 1. 



(63) 



From the results of the previous section, we know that 
bk = c\(hj J) K a'l where ci is a constant (equal to 
(co/8(d — l)) 2 for the lattice case). Now n is greater 
than zero for both the lattice and the continuum. Thus, 
for small h, (|6^) becomes b k = 1 so the final level of 
decorations, k, is determined by 



a k = L/n k = 
This gives the cutoff k as 



lni+ ilnci/i K 
Inn 



(64) 



(65) 



This result should be approximately valid for large L. 

Finally, all the changes to the free energy from all 
the various decorations are summed to arrive at a to- 
tal decoration. At stage i we have decimated the length 



6. Scaling and criticality in d — 2 and d = 2 + e 

The preceeding results can now be used to construct 
the RG transformation in d = 2 + e. This allows us 
to investigate the marginality breaking effect of domain 
wall roughening in two dimensions, and to investigate 
the behaviour of the h = fixed point as we move away 
from d = 2. So we obtain the two dimensional critical 
behaviour and also an idea of the critical properties in 
higher dimensionalities. 

Using (^36|) and the free energy forms found in the pre- 
vious sections, a full decoration of the domain walls may 
be done in general dimension d > 2. For the sake of 
generality we wi ll us e the general form of the wall free 
energy given by ( 32a ) , though with restored generic con- 
stants as well as exponents, which can be read off from 
earlier results for specific cases. Thus, the general form 
of the wall free energy is 



L e^ n i(d-i- 



(67) 



i=0 



From the values of 9 given in Table IV, the summand is 



i(d-i-e) _ 



1 + iielnn , continuum 



4 i(2d-4) = l + - 2elnj 



lattice 



for d = 2 + e. To further simplify matters we define the 
coefficient 



continuum 
lattice 



(68) 



Then, the total free energy of the wall decorations may 
be concisely written as 



6? 



■I'ot 



,|'^| L" 



(k+l) + e<f) (Inn) 



k(k + l) 



(69) 

This result may be compared both with that of Villain^ 
and Grinstein and Mala who find 5T Tot - h^LhiL for 
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d = 2 using a continuum interface model, and with the 
result of BinderEl who finds &F Tot - h 2 L\nL in a 2d 
lattice calculation. Both these results agree with j69| ) 
with the appropriate choice of the exponents 7 and 9 
(see Table 0). 

The wall decorations have the effect of reducing the 
effective surface tension of the domains. As a result, the 
average domain size will decrease. We can quantify this 
result by rewriting the free energy of the entire strip as 



N 



T={ — ) (2JL d - 1 - c hy/E L L d - 1 + 



Tot 



(70) 



Here, we used the fiat wall result (||) as a model, but 
noted the fact that there is, on average, a further contri- 
bution of ^JF Tot per domain wall from the decorations. 

Extremalizing this free energy with respect to gives 
the average domain size as 



cp h 
J 



-d-l 



5T 



Tot 



JL 



d-1 



(71) 



As expected, one sees that the domain size is reduced by 
the effect of interface roughening. 

This correction to the correlation length also has 
implications for the phcnomcnological RG equation. This 
equation obtained from (|J) with the use of (|7l|) and (^|) is 
now quite complex, but it can be simplified by linearizing 
in e and using small h/ J. After some tedious algebra, one 
finds the following RG equation: 



Inn 2 



(72) 



and we have let J = 1 for simplicity. Be- 
sides the trivial fixed point at h = 0, ( |72|) has an unstable 
fixed point at 



elnn 
2c 



l/ 7 



2n1/ 7 



(73) 



which remains in the small h region of validity of our 
description since e is regarded as small. The eigenvalue 
of (173) is given by 



Bti 
~dh 



1 



= 1 + - 7 eln6 + 0(e 2 ) = b 1 '". 



h=h' 



This yields the critical exponent 

2 

v = — 



(74) 



which characterises the scaling of the bulk correlation 
length, £ ~ \h — h*\~ u . Notice that the arbitrary scale 
variable n disappears from the critical exponent v, which 
is a universal quantity, though it remains through a weak 
logarithmic dependence in h* and other non-universal 
quantities. 



If we attempt to let e — > to extrapolate back to the 
h = fixed point, we see that v diverges. This is an indi- 
cation that the divergence of the correlation length of the 
infinite two dimensional system is faster than any power 
law. This divergence can be obtained from equation ( |7^ ) 
at e = 0: 



h! = h 1 + ch 



Mb 
Inn 

This is equivalent (at small h) to 
/ Inn \ 1 



exp 



\ic(h') 



( Inn \ 
b eXP {^) 



(75) 



(76) 



This is the same form as the standard rescaling of a cor- 
relation length £(h') = 6 _1 £(/i), so the correlation length 
is 



(77) 



close to the fixed point h* = 0, where A is the non- 
universal constant Inn/ (07). 

The results of these calculations are summarized in 



Table III 



Here the exponent is 7 = 4/3 for the continuum or 7 = 
2 for the lattice. r-^The exponential growth of £ was also 
found by BinderEl His calculation was done explicitly 
on a 2d lattice and he found a form for £ identical to 
ours with 7 = 2. 



7. Conclusions 

This concludes the T = investigation of the effects 
of fluctuations in the domain wall shape. The analysis of 
the RFIM on a lattice and in a continuum has led to a dif- 
ference in the critical behaviour of these two cases which 
arises from the difference of their wandering exponents. 
This has been seen previously in other works referred to 
above. One may speculate that this difference is due to 
the different effect of surface tension in the lattice versus 
the continuum. Interfacial roughening effectively acts to 
reduce the surface tension in the domains and is the cru- 
cial factor in breaking the marginality of h in d = 2. This 
surface tension is coupled to the d — 1 dimensional area 
of the domain wall and this is clearly different between 
lattice and continuum cases because of the stepped char- 
acter necessary on a lattice. Exactly how this mechanism 
comes into play in the RFIM is not clear, however. 

The procedure was to perform decorations on all length 
scales to find the effect of domain wall roughening on the 
free energy and domain size. These determine the phe- 
nomenological renormalization group equation and hence 
the behaviour near the lower critical dimension. This re- 
vealed the evolution of the fixed point as it moves from 
h = to h = h c ^ for d > 2. 
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As with the flat wall analysis, both the resulting scal- 
ing and the ingredients which lead to it (size and field 
dependences of free energy and domain size, and now 
also roughening) can be directly che cked by numerical 
work, and this will be done in section III C . 



In the next section, we extend our results to include 
low temperature thermal effects. 



B. Low Temperature Theory 

The effect of the zero temperature wall roughening on 
the surface tension is (by (|70|)) equivalent to the replace- 
ment of J by 



J=J + 



ST 



Tot 



(78) 



Cons eque ntly, the lo w temperature analyses given in sec- 
tions [I B and [I C 2 have to be modified by the replace- 

D-(F 



ment J -> J, in equations (|20|)-([24|) and 

The roughening also affects the entropy since it in- 
creases the number of possible configurations. The en- 
tropy change can be obtained crudely by a slight gener- 
alisation of the counting used to derive flog) , or by a more 
complete analysis given el sewh erec2l. 

As remarked in section II B , the zeroth order thermal 
scaling (^3|) is not marginal in d > 1, and so the modifi- 
cations just discussed do not change the RG thermal flow 
directions for h s mall an d d near 2. So the e-expansion 
results of section [II A 6) can be combined with those of 
section [I B to infer the flow diagrams shown in Figure [l4| 
for d = 1, 2, 2 + e. For d = 2 + e this implies a second 
order phase transition at the phase boundary shown, in- 
side which (i.e. at small h, T) long range order occurs for 
h =/= 0. This is co nsist ent with previous renormalization 
group discussionsc2l,E2l. 

The extended theory outlined in the first two para- 
graphs above provides the low temperature free energy 
modifications produced by rough walls. 

This suggests the possibility of distinguishing wall 
roughening effects by comparing numerical data on free 
energy with this theory extended as just outlined. This 
comparison has been made with the transfer matrix free 
energy data discussed in section II C 2. Despite the excel- 



lent quality of the data, and its remarkably good agree- 
ment with the theory in zeroth order (flat walls), it was 
not found possible to extract wall roughening effects from 
it. This is because the statistical fluctuations in the data, 
although extremely small, are nevertheless greater than 
the roughening contributions in the ranges of parame- 
ters appropriate for the validity of the theory (ensuring 
£ S> L). So, we turn to other more sensitive comparisons. 



C. Numerics - Correlation and roughening 
characteristics from the Max-flow ground state 
algorithm 

This section describes the extraction of domain wall 
roughening characteristics from data obtained using the 
max-flow algorithm for constructing ground states. 

One comparison uses domain size correlation length 
data, which is much more sensitive than the free energy 
to domain wall roughening, particularly in wide strips 
because of the crucial role of roughening in the marginal 
dimension d — 2. The other involves direct measures 
of the wall wandering. These measures are the average 
over many walls, (for a given h and base scale L) of the 
maximum height (wandering excursion) b max , or of the 
root mean square height b rms . 

Both binary and Gaussian random field distributions 
were used. As expected, the Gaussian distribution gives 
smoother dependences, and has provided the most useful 
data. But, the binary results are interesting. Figure [l5| 
gives numerical data for the configurationally averaged 
domain size. A "stepping" of the data as a function of 
h is apparent, particularly for smaller L. This is pre- 
dicted by the flat wall description (section II A 5|) equa- 
tions @, (0) since in (0), q = (2 + [(2 JL) /h\ ) _1 
involves the integer part function |_ J . There is also a cer- 
tain discreteness implied in the use of the transfer matrix. 
For comparison, the flat wall theory is shown in the Fig- 
ure. The flat wall random walk analysis also compares 
very satisfactorally with corresponding domain size data 
(Figure [To]) from the Gaussian distribution. Here, the 
data contain no steps. If one repeats the random walk 
analysis with a Gaussian random field distribution, the 
"integer part" function does not occur, and the trans- 
fer matrix becomes a convolution operator. Since the 
eigenfunctions are again plane waves, the final result for 
the correlation length is the same, except that it has no 
"integer part" function, and thus no "stepping". But, 
though the Gaussian distribution must therefore set the 
large length scale behaviour, in energy (Figure ||), do- 
main size (above) and probably in time scales for relaxing 
on given characteristic scales ("dynamics"), there can be 
real differences for these quantities between binary and 
Gaussian cases at smaller scales. 

The domain size data can be fitted to the flat wall 
theor y or to the ana lysis including wall roughening (sec- 
tions III A 5| , III A 6 ). Equation (Q) gives the basic flat 
wall description, in which for strips [d = 2), H/, ~ Lh~ 2 , 
and one fitting coefficient (cq) occurs. The fuller ran- 



dom walk version of the flat wall theory (section II A 5) 



is equivalent, for small h/L, but more generally it gives a 
El different by a "discreteness factor" (l + 0(h/L)). The 
wall roughening analysis instead corrects the flat wall 5 l 
by a factor (l + 0(h~< L - d+1 lni)). This should be amal- 
gamated with the discreteness factor evident from the 
RW analysis. Then the domain wall roughening theory 
gives the best quality of fit to the max-flow data on Si. 
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For instance, considering a Gaussian random field distri- 
bution, Table shows the decreased statistical error x 2 
associated with the addition of this domain roughening 
term. Even though this extra term provides an extra 



free constant to the fit (c 



it does not necessarily 



decrease x 2 ■ Indeed, it worsens the fit to the zero'th or- 
der theory since the higher order "discreteness term" has 
been neglected. When this term is ignored, the roughen- 
ing only has a correlation coefficient of r — 0.218. How- 
ever, when this discreteness term has been included from 
the more complete random walk analysis, the wall rough- 
ening is seen to decrease the value of \ 2 . It then has a 
much more significant correlation coefficient of r = 0.816. 

The max-fiow data in the wall decoration variables 
brms, b m ax give much more conclusive evidence for the 
roughening effects and provide quantitative estimates of 
exponents. Data for the /i-dependence of b rms for L = 8 
are shown in Figure [l9] for both binary and Gaussian 
distributions. Again, the jagged character of the binary 
result is evident, even at this moderately large L. So the 
remaining discussion is confined to the Gaussian case. 
To compare with the prediction ( |62b| ), a log-log plot of 
the Gaussian b rms versus h is given for L = 2, . . . , 11 in 
Figure [l8| A very similar plot is obtained also for b max: 
though the absolute size of b ma x is larger, by up to a 
factor of 4. The departure from linearity of the log-log 
plot at small h is almost certainly a lattice effect. That 
requires b>l, and using e.g. (|38| ) ( with Co = 1.8) to 
estimate where b ~ 1 suggests that ( |62b| ) should break 
down for h/J < 5L -1 / 2 , which is in qualitative agreement 
with what is se en. Similarly, a breakdown of the simple 
theory giving ( 62b ) is expected when b (strictly b rnax ) 
becomes comparable to Hj,, i.e. (using ( p8| ) and (^)) 
for h/J>1.5. This L-independent cutoff is consistent 
with the departures from linearity seen at the largest 
fields in the log-log plot. The scaling window between 
the upper and lower h departures is quite wide for the 
largest L's, and consequently they should give the mos t 
accurate field exponent n from power law fits to (321) 
within the scaling window. The results, from both b rms 
and b max are shown in Figure [l^ and are consistent with 
k w 2.1±0.3. Figure p"7| shows the L-dependent prefactor 
of h K in b rms and b max . In each case, the prefactor is 
very close to linear in L i n th e larger L scaling regime, 
consistent with £ = 1 in (62h). The behaviours seen in 
Figure |l8|, and in the corresponding plot for b max , as well 
as the associated exponents k, ( are therefore consistent 



with the analysis in section III A 



IV. CONCLUSIONS 



is removed by wall roughening effects (Section III). The 
resulting phcnomcnological renormalization group trans- 
formation has been used to obtain the critical properties. 
In particular, the special critical behaviour of the two- 
dimensional correlation length and the ph ase diagram in 
d — 2 + e have been given (Section [IIIA6; ). 

The wall roughening ingredient has its own scaling 
characteristics (the exponents 7, 9, k, £) which have been 
here discussed usingjJpoth a simple analytic method 
(following NattermanE 2 ]) and a field theoretic approach. 
These exponents have a direct bearing upon the RG 
transformation and hence on the scaling properties. 

An essential element to this study has been the support 
that numerical studies have provided for the analytic de- 
scription. The comparison between numerical data and 
analytic predictions was made in the strip geometry so 
that the critical ingredients of the theory could be di- 
rectly investigated. These ingredients include the domain 
size, the free energy, and the roughening characteristics, 
as functions of h, T and strip width L. 

The numerical approaches employed were transfer ma- 
trix and Monte Carlo techniques, and a ground state 
(max-flow) algorithm. The first two give free energy 
and domain size data confirming the basic flat wall pic- 
ture through: (i) the h 2 /L dependence of the T = flat 
wall energy at low h, and its low temperature generalisa- 
tion; (ii) the domain size distribution and the correlation 
length-domain size relationship. 

The numerical ground state data, obtained via the 
max-flow algorithm, directly verify the basic validity of 
the domain size predicted by the (flat wall) random walk 
analysis, but they also provide evidence for the roughen- 
ing effects in the domain size. The max-flow data on the 
wall decoration variables b rms , b max are the most con- 
clusive evidence for roughening effects. They give power 
law scaling in a window of the size predicted by the the- 
ory, sufficient at large L to give quite accurate values 
for roughening exponents, in good agreement with the 
(lattice based) theory. 

The conclusion is that the theory contains the correct 
ingredients, and the numerical data provide quantitative 
confirmations of the way they enter into the theory. 

We have not made corresponding numerical investi- 
gations of the d — 3 case, where it is of course more 
difficult to obtain data for the range of L's we needed 
here. However, we hope the present work may stimu- 
late such an effort, which (e.g. using the ground state 
algorithml-fpnld sort out the so far unresolved critical 
behaviounljilil. Another extension, presently under con- 
sideration, is to the kinetic behaviour, for which the free 
energy scaling discussed here is an essential ingredient. 



In this paper we have developed a domain scaling de- 

scription for the random field Ising model by exploiting 
a bar geometry. As in earlier work of Villairo, Grin- 
stein and Mai and Bindero, the marginality of the basic 
flat wall picture at the lower critical dimension di = 2 
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FIG. 1. The bar geometry, showing the flat domain walls 
that separate regions of up and down spins. The average 
domain size El is the average of the Hi's. 



FIG. 2. Illustration of how the random field configuration 
may be viewed as a random walk in 1 + 1 dimensions. We use 
L — 3 here. 



FIG. 3. Log-log plots of the excess free energy at T = 0.1J. 



FIG. 4. A comparison of the free energy as determined by 
the transfer matrix data, and that determined by our low 
temperature theory using the fitted constant Co = 3.406, at 
T = 0.1J and T = 0.3J. 



FIG. 5. The free energy as a function of h for T/J = 0.1 
and L € (2, 7). These data were determined by the transfer 
matrix method, and illustrate the three regimes of h depen- 
dence 



FIG. 6. The average energy of the L = 1 and L = 4 
RFIM, as determined by both the Monte Carlo algorithm 
(data points) and the zeroth order theory (lines). The sym- 
bol O corresponds to h = 0, — to h = 0.1J, □ to h = 0.2J 
and x to ft = 0.5J. 



FIG. 7. The correlation length £l plotted against the do- 
main size El for L — 2. We also include the predicted relation 
between these quantities for comparison. Plots for L = 1,4 
are essentially the same. 



FIG. 8. The distribution of domain sizes according the the 
measure Sl as found by the Monte Carlo routine when L = 1. 
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FIG. 9. The domain size E L plotted as a function of tem- 
perature for L = 1 and L — 4. The data points represent 
the MC data whereas the lines show the result of solving the 
flat wall theory. The symbol O corresponds to ft = 0, — to 
ft = 0.05J, □ to ft = 0.1J, x to ft = 0.2J, and A to ft = 0.5J. 
The flattening of the data at low T is due to the finite length 
of the simulated system. 

FIG. 10. The typical domain size at T = from the 
max-flow algorithm. The points represent the data from the 
max-flow algorithm, and the solid lines represent the theoret- 
ical form with a fitted multiplicative constant. 

FIG. 11. First order decorations to a domain wall on a 
lattice (left) and in a continuum (right). Here, d = 3. The 
decorations have height b and width a. 

FIG. 12. Plot of the free energy exponent 6 and the wan- 
dering exponent £ for the wall decorations on a lattice and in 
a continuum. 

FIG. 13. The process of placing decoration upon decora- 
tions, such that decorations occur on all length scales. This 
process has a cut-off length scale when the size of the deco- 
rations become on order of the lattice constant. 

FIG. 14. The RG flow diagram for the RFIM as predicted 
by our finite size scaling analysis. 



FIG. 19. The rms domain wall width as a function of ft/ J 
for L = 8. This plot emphasizes the difference between the 
binary and Gaussian random field distributions, showing the 
stepped nature of the former as more excitations become 
favourable at lower h/J. 



TABLE I. Fits to the free energy data from the transfer 
matrix 



L Slope 


Error 


Q 


2 1.792 


0.007 


0.0004 


3 1.848 


0.009 


0.33 


4 1.857 


0.009 


0.79 


5 1.885 


0.011 


0.96 


6 1.901 


0.012 


0.95 


7 1.904 


0.012 


1.00 


8 1.945 


0.015 


1.00 


9 1.958 


0.014 


1.00 


TABLE II. Fits of the free energy to the transfer matrix 


data at finite temperature 






L N 


2 

co X 


Q 


2-9 120 


1.774 457 


1( j-41 


2-7 90 


1.761 375 


io- 37 


2-5 60 


1.745 288 


1Q -32 


4-9 90 


1.846 81 


0.686 



FIG. 15. The domain size as a function of h/J at T = 
determined by the max-flow algorithm (data points) , and the 
fiat wall analysis (lines) for a binary field distribution. 

FIG. 16. The wall roughening exponent k as determined 
from power law fits to b rms and b max . For large L, these tend 
towards 2. 



TABLE III. Main results of the Flat Wall analysis 



Dimensions Fixed Point Bulk Correlation Length 

2 ft^O 1 - exp(A/ft 7 ) 

2 + e ft* = (^) 1/7 £~(ft-ft*)W 



FIG. 17. The wall roughening term or as determined from 
power law fits to 6 rms and b max . For large L, these tend 
towards a linear dependence on L, consistent with £ = 1, as 
shown by the lines of best fit. 



FIG. 18. The root-mean-square domain wall width as a 
function of h/J for various length scales L. The use of a 
log-log plot shows the scaling window where the plots are 
linear. 
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TABLE IV. Comparison of Approaches to Wall Decoration Problem 







d < 2 




a ^ z 






Lattice 


Simple 


Field 


Lattice Simple 


Field 






Continuum Theory 


Continuum 


Theory 


c 


3-d 


3-d 


N.A. 


3-d (5-d)/3 


(5 - d)/3 


e 


3-d 


1 


N.A. 


3-d (d+l)/3 


(d+l)/3 


7 


2 


2 


N.A. 


2 4/3 


4/3 


K 


2 


2 


N.A. 


2 2/3 


2/3 






TABLE V. 


Results of fits to the T 


= domain size data 








Flat Domain Walls 


Roughened Domain Walls 




Zeroth Order 


Random Walk 


Zeroth Order 


Random Walk 


x 2 


4.215 




0.607 


4.447 


0.404 
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